Model independent determination of the dark matter mass from direct detection 

experiments 
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Determining the dark matter (DM) mass is of paramount importance for understanding dark 
matter. We present a novel parametrization of the DM speed distribution which will allow the DM 
mass to be accurately measured using data from Weakly Interacting Massive Particle (WIMP) direct 
detection experiments. Specifically, we parametrize the natural logarithm of the speed distribution 
as a polynomial in the speed v. We demonstrate, using mock data from upcoming experiments, that 
by fitting the WIMP mass and interaction cross-section, along with the polynomial coefficients, we 
can accurately reconstruct both the WIMP mass and speed distribution. This new parameterisation 
represents a significant improvement over previous methods as it provides an accurate, unbiased 
reconstruction of the WIMP mass even for non-trivial underlying speed distributions, demonstrating 
that the WIMP properties can be reconstructed from direct detection experiments without prior 
knowledge of the distribution function. 

PACS numbers: 07.05.Kf,14.80.-j,95.35.+d,98.62.Gq 



Dark matter (DM) has thus far only been detected 
through its gravitational interaction with Standard 
Model particles and its particle nature is not yet known. 
Weakly Interacting Massive Particles (WIMPs) are a 
good DM candidate as they are generically produced in 
the early Universe with the required abundance and Su- 
persymmetry (SUSY) provides a concrete well-motivated 
WIMP candidate in the form of the lightest neutralino 
(e.g. Ref. [l|). WIMPs can be detected directly in the lab 
or indirectly via their annihilation products. They can 
also be produced at particle colliders such as the LHC. 
Direct detection experiments [2|, [3[ aim to observe nuclear 
recoils produced by WIMPs passing through terrestrial 
detectors (such as XENON100 [4] or CDMS @). While 
WIMPs are within the reach of near future direct detec- 
tion experiments their convincing detection is likely to 
require consistent signals (i.e. with the same WIMP pa- 
rameters) from direct, indirect and collider searches (e.g. 
Ref. [6]). An accurate determination of the WIMP pa- 
rameters would also be instrumental in constraining the 
parameter space of SUSY models (e.g. Ref. [7[). 

Measuring the rate of nuclear recoils as a function of re- 
coil energy should in principle allow the WIMP mass and 
interaction cross-section with nucleons to be extracted. 
However, the analysis of direct detection data requires 
assumptions to be made about the astrophysical distri- 
bution of the DM within the Milky Way halo. The local 
velocity distribution, /(v), encodes the speeds of DM 
particles and determines the recoil energies observed in 
experiments. Direct detection experiments usually as- 
sume the simplest possible model for the Milky Way halo 
(referred to as the Standard Halo Model). This model 
assumes that the halo is isotropic, has a density distribu- 
tion p(r) oc r -2 and has reached a steady state, in which 
case it has a Maxwell-Boltzmann velocity distribution in 
the Galactic frame. However, high resolution N-body 



simulations [8|, [9[ suggest that the true distribution func- 
tion may be non-Maxwellian. In particular there may be 
features in the high-speed tail of the speed distribution 
from particles that are not yet virialised [Toj |. Further- 
more, the effect of baryons on the DM halo is not yet 



fully understood. For example, some simulations [111 Il2] 
show evidence for a dark disk (DD) which corotates with 
the baryonic disk and which may significantly affect the 
direct detection rate fl3j . 

These uncertainties in the velocity distribution can 
lead to an order of magnitude uncertainty in estimates 
of the interaction cross-section [14J and significant bias 
in the recovery of the WIMP mass from direct detec- 
tion data [TBI. It is therefore imperative to account for 
this uncertainty in the current and future analysis of 
such experiments. Several different approaches have been 
proposed. One option is to simultaneously fit Galac- 
tic parameters, for example the lag speed (the speed 
of the Solar System with respect to the peak of the 
WIMP speed distribution), the velocity dispersion and 
anisotropy, alongside the WIMP mass and cross-section 
0, [TBI, Eg • However this method assumes that the Milky 
Way halo is fully equilibrated and requires that it can be 
described by one of a relatively small class of models. Pe- 
ter proposed a more model-independent approach where 
f(v) is parameterized as a series of constant bins in veloc- 
ity space [15]. Subsequent studies have explored param- 
eterizing the momentum distribution 18], while Ref. [~L9| 



parameterises f(v) in terms of integrals of motion. How- 
ever, these approaches still have significant shortcomings 
and either result in a substantial bias in WIMP param- 
eters [TB[_ 18] or assume that the WIMP mass is already 
known]l9t 

In this letter we present a new model-independent 
parametrization method for f(v) which allows the ac- 
curate, unbiased reconstruction of the WIMP mass. In 
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order to demonstrate the robustness of the method, we 
generate mock data sets for three proposed experiments 
(modelled on detectors which are in development) using 
several benchmark speed distributions. 

The event rate R per unit nuclear recoil energy Eji in 
a DM detector can be written as [fll 
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where po is the local DM mass density, m x is the WIMP 
mass, a p is the WIMP-proton cross-section (at zero- 
momentum transfer) and \i xv is the WIMP-proton re- 
duced mass, /i XP = (m x m p ) / (m x + m p ). Here, we have 
considered only the spin-independent contribution to the 
rate, which dominates for targets with large atomic mass 
number, A. Assuming that the coupling to protons and 
neutrons is identical, this leads to an A 2 enhancement in 
the scattering rate Finally, the loss of coherence due 
to the finite size of the nucleus is accounted for by the 
form factor, F 2 (Er), which we take to have the Helm 
form 

The local velocity distribution of Galactic WIMPs en- 
ters into the event rate through r](v m i n ), the mean inverse 
speed: 
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where /(v) is the normalised velocity distribution in the 
Earth frame and v m i n is the minimum WIMP speed re- 
quired to produce a recoil of energy Er. For a detec- 
tor of nuclear mass, tun, and reduced WIMP-nucleon 
mass, fi X N = (m x rriN)/(m x + rajv), this minimum re- 
quired speed is 
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Detector Target Xe 


Ge 


Ar 


Nuclear Mass, m^/amu 131 


73 


40 


Detector Mass, radet/kg 1000 


100 


1000 


Exposure Time, t exp /days 152 


274 


912 


Energy Range, [2,30] 


[10,100] 


[30,130] 


[E min , #max] /keV 
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— - SHM + Dark disk 



TABLE I. Parameter values for the three mock experiments 
used in this work. The meanings of the experimental param- 
eters are described in the text. 
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FIG. 1. The speed distribution for the three benchmark mod- 
els defined in the text: SHM (solid blue), SHM+DD (dot- 
dashed green) and stream (dashed red). 
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In order to generate mock data sets, we use three hy- 
pothetical experiments which are representative of those 
currently in development (namely XENON1T [21], Su- 
perCDMS [lH and WArP [23]) and we assume for con- 
creteness that our proposed detectors have perfect energy 
resolution and negligible backgrounds. Each experiment 
is then characterised by the mass of the target nucleus 
in 7v i the total detector mass radet? the effective exposure 
time £ e xp and energy window accessible to the experi- 
ment [E'min? -E'max]- Information about detector efficiency 
is accounted for in the value of t exp . The values used for 
these parameters in the current work are shown in Table 
ID These have been chosen to closely match those used in 
Ref. [15] and [18], though with increased exposure times. 

We consider three possible underlying distribution 
functions. The first is the Standard Halo Model (SHM), 



which we model as a Gaussian in the Galactic frame with 
width a v = 156 km s _1 and lag speed ^i ag = 230 km s _1 . 
The remaining two extreme distributions lead to signif- 
icantly different results to the SHM and one would ex- 
pect them to be difficult to reconstruct. These are: a 
DM stream (v\ ag = 400 km s _1 , a v = 10 km s _1 ) and 
an SHM with a 50% admixture of a dark disk (v\ ag = 
a v = 50 km s _1 ). These are referred to as 'stream' 
and 'SHM+DD' respectively. All three distributions are 
shown in Fig. [TJ 

The experiments used here are typically insensitive to 
low mass WIMPs (m x < 10 GeV), while higher mass 
WIMPs (m x > 100 GeV) lead to a degeneracy in the 
m x — <j p plane which has been explored thoroughly in the 
literature (see Ref. [TH, H^, for example). We therefore 
consider a single benchmark WIMP mass of 50 GeV in 
the current work. 

We use a fiducial cross-section of a p = 2 x 10 -45 cm 2 , 
which is at the current bound from XENON100 [4]- How- 
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ever, this bound was calculated assuming a Standard 
Halo Model and, in any case, the value of a p is de- 
generate with the local DM density po? which is typi- 
call y p resented with substantial uncertainties (e.g. Ref. 
25|-l2S[)- We therefore choose to fix the cross-section at 
the above fiducial value and fix the local DM density at 
po = 0.3 GeV cm" 3 . 

For each set of benchmark parameters we generate 
a single realisation of data. The impact of Poissonian 
statistics on reconstructions has previously been studied 
in detail in Ref. [29|. Using these benchmark values, the 
total number of events observed across all three experi- 
ments ranges from 176 in the SHM+DD case to 386 for 
the stream. 

We assume no prior knowledge of the speed distribu- 
tion and attempt to construct a parametrization which 
will allow us to fit the mock data sets in a model indepen- 
dent way. The one (rather weak) assumption we make 
about the speed distribution is that f(v> v max ) = 
for v max = 1000 km s _1 . The choice of the value v max is 
somewhat arbitrary, but this particular value is conserva- 
tive as the Galactic escape speed is significantly smaller 
than this [8|. 

We note that, by definition, the speed distribution 
is everywhere positive f(v) > 0, motivating us to 
parametrize the logarithm of f(v) as a polynomial of de- 
gree N in v. We therefore write 



N 
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meaning that the full directionally-averaged speed distri- 
bution is 



fi(v) = v 2 exp - 



N , 
^CLkPk ( 

k=o ^ max 



subject to the normalisation condition 

fi(v)dv = 1 . 



(5) 



(6) 



We have chosen a basis of shifted Legendre polynomials 
Pk of degree k = {0, 1, TV} for the parametrization. In 
theory, any basis may be used, but in practice choosing 
polynomials which are orthogonal over some finite range 
of v improves the behaviour of the coefficients a&. By 
varying TV, the number of terms in the polynomial, we 
can accommodate features in the distribution function, 
such as multiple components. By widening the priors on 
the values a/~, we can accommodate sharper individual 
structures. 

We explore the posterior distribution for the param- 
eter space |m v , cr p , ja^ II using the publicly available 
MultiNest package [5o|, [H[ as a generic nested sam- 
pler. We use an unbinned likelihood function (as intro- 
duced in Ref. [32[) and use the following sampling pa- 
rameters: N live = 10,000, eff = 0.3, tol = 1.0e-4. 
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FIG. 2. 68% and 95% confidence contours for the WIMP 
mass m x and cross-section a p obtained using the speed 
parametrization described in the text. Results are shown 
for three different underlying speed distributions: SHM (solid 
blue), SHM+DD (dot-dashed green) and stream (dashed red). 
The true parameter values are shown as a black cross. 



We use log- flat priors on m x G [1, 1000] GeV and a p G 
[1,1000] x 10 -46 cm 2 , linearly-flat priors on \ai\ < 50 
and go up to degree N = 4 in the basis polynomials. For 
the stream benchmark, we allow \ai\ < 500 and for the 
SHM+DD distribution, we extend the basis up to N = 9. 
We can straightforwardly check the robustness of a given 
reconstruction by increasing N and widening the priors 
on a^- 
In Fig. [2j we show the 68% and 95% confidence con- 
tours for the reconstruction of the mass and cross-section, 
obtained from the 2-dimensional profile likelihood. In 
the case of the SHM and stream, the WIMP mass and 
cross-section are both well recovered, lying within the 
68% contour. In the case of the SHM+DD distribution 
function, however, there is a significant bias towards a 
smaller cross-section, with a discrepancy of around 30% 
between the underlying and reconstructed cross-section. 

Table [HI shows the best fit WIMP masses and error es- 
timates derived from these profile likelihoods. In contrast 
to the bias observed in the reconstruction of the cross- 
section, the WIMP mass is accurately reconstructed in 
all three cases, with the true mass lying within the l-cr er- 
rors. For comparison, Table HT1 also shows the value of the 
mass reconstructed using the speed binning method of 
Peter [15] and the momentum binning method of Ref. [HI 
for the stream benchmark. Both of these methods show 
a bias towards lower masses and significantly underes- 
timated errors. The new parametrization method pre- 
sented here is seen to perform significantly better in re- 
covering the WIMP mass. 



Speed distribution Reconstructed mass (GeV) 
Standard Halo Model (SHM) 50.0^ 2 2 6 
SHM + Dark disk 44.2±f;£ 
Stream 44.7t^ 
Stream (using 5 speed bins) 29.3^5'q 
Stream (using 5 momentum bins) 38.2lJ>'.3 



TABLE II. Reconstructed mass calculated from the profile 
likelihood, using the parametrization presented in this letter. 
Also shown is the mass reconstructed from the stream distri- 
bution using the speed binning [15|] and momentum binning 
[HI ] methods for comparison. 



In Fig. [3j we show the reconstructed speed distribu- 
tions for the three benchmarks. The 68% credible regions 
(shaded pink) are obtained by marginalising over the val- 
ues of f(v) - derived from the posterior distribution of the 
{dk} - at each value of v. The errors in f(v) across dif- 
ferent speeds are therefore strongly correlated. However, 
such reconstructions should be representative of the un- 
derlying form of the distribution (solid line). Also shown 
is the best fit form of f(v) in each case (dashed line). 

For all three benchmarks, the speed distribution is 
largely well recovered. The SHM benchmark recovers 
a broad distribution, while the 50% dark disk produces 
a distinctive double peak structure. The stream bench- 
mark results in a sharp localised peak within ~ 20 km s _1 
of the true peak. However, the reconstruction is typically 
poorer for low and high speeds. This is particularly not- 
icable in the SHM+DD benchmark for v < 100 km s _1 . 
This is due to the fact that for a WIMP of mass 50 GeV 
and the energy thresholds given in Table HJ the small- 
est speed to which the detectors are sensitive is v ~ 
92 km s _1 . Below this speed, the distribution function 
is poorly constrained. Similarly, at high v, the upper en- 
ergy limit E max , combined with form factor suppression 
of the event rate, leads to low sensitivity. 

Due to the energy thresholds of the experiment, we 
cannot constrain f(v) across its entire range and there- 
fore we do not know what fraction of DM particles lie 
within reach of the experiments. Generically then, it is 
not possible to reconstruct cr p , and the bias in <j p shown 
in Fig. [2] is a problem inherent to all model independent 
approaches. In the case of the SHM+DD benchmark, 
only 72% of WIMPs lie within the sensitivity range of 
the experiments, leading to a reduction of around 30% 
in the effective cross-section. In the case of the SHM 
and stream benchmarks, almost all WIMPs are accessi- 
ble to the experiments, leading to an accurate reconstruc- 
tion. However, we do not know a priori if this is the case 
and therefore it is only possible to obtain a lower bound 
on the cross-section. Nonetheless, this ignorance of the 
cross-section does not preclude us from accurately recon- 
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FIG. 3. Reconstructed speed distributions for the SHM (top), 
SHM+DD (middle) and stream (bottom) benchmarks. The 
true distribution is shown as a solid curve, while the best fit 
reconstructed value is shown as a dashed curve. Also shown 
are the marginalised 68% credible intervals (shaded pink). 



structing the DM mass, as we have seen. 

In summary, accurately measuring the WIMP mass 
from direct detection data is crucial in order to convinc- 
ingly detect, and characterise, DM. We have presented a 
new parametrization of the DM speed distribution which 
does so in a model independent way. While this leads 
to an unavoidable uncertainty in the cross-section, it 
does allow the speed distributions to be well constrained 
within the range of sensitivity of the experiments. Most 
importantly, however, the new parametrization is able to 
robustly and accurately recover the WIMP mass, even for 
non-trivial underlying speed distributions, representing a 
significant improvement over previous methods. 
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